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ABSTRACT 

We show how the morphological analysis of the maps of the secondary CMB 
anisotropics can detect an extended period of "smoldering" reionization, during which 
the universe remains partially ionized. Neither radio observations of the redshiftcd 
21cm line nor IR observations of the redshifted Lyman-alpha forest will be able to de- 
tect such a period. The most sensitive to this kind of non-gaussianity parameters are 
the number of regions in the excursion set, N c i, the perimeter of the excursion set, P g , 
and the genus (i.e. '1 - number of holes') of the largest (by area) region. For example, if 
the universe reionized fully at z — 6, but maintained about 1/3 mean ionized fraction 
since z = 20, then a 2 arcmin map with 500 2 pixel resolution and a signal-to-noise 
ratio S/N = 1/2 allows to detect the non-gaussianity due to reionization with better 
than 99% confidence level. 

Key words: cosmic microwave background - cosmology: theory - cosmology: large- 
scale structure of universe - galaxies: formation - galaxies: intergalactic medium 



1 INTRODUCTION 
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the physical conditions in the intergalactic medium (IGM) 
shortly after the epoch of reionization. Similar future obser- 
vations will increase the amount of observational data multi- 
fold and will provide critical constraints on the theories of 
galaxy formation. 

However, if we want to go beyond the epoch of reion- 
ization and to study the pre-overlap stage during which ion- 
ized Hll regions expand into still neutral low density gas, 
we need to use different bands of electromagnetic spectrum. 
Next Generation Space Telescope (NGST) will provide valu- 
able clues on the earliest episodes of galaxy formation from 



Pen 2001; Gnedin & Jaffe 2001). While such a measurement 



is also in the future, observations of the secondary CMB 
anisotropies can provide constraints on the physical condi- 
tions in the IGM that are not accessible by other means. 
In fact, because CMB anisotropies are sensitive to the total 
Thompson optical depth, they can probe low levels of ion- 
ization in the mostly neutral gas - which cannot be done by 
radio observations of the redshifted 21 cm line. 

Several scenarios have been proposed recently in which 
the universe undergoes a protracted episode of incomplete 
ionization prior to full reionization, cither due to early 
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infra-re 1 observations, although a relationship between the 



first galaxies and the properties of the IGM at high redshilts 
will not necessarily be easy to determine from such obser- 
vations. Radio observations of the redshifted 21 cm line of 
neutral hydrogen might be capable of measuring the pre- 
reionizati on signal ( MadMi^M^^^n^^^^es^99^Jro^^ei 
al. 2000| ; pev et al. 2002j ; jCarilli. Gnedin, &r Owen 2002[ ), 
although such a measurement will be at the very edge of the 
capability of the next generation radio instruments such as 
Low Frequency Array (LOFAR). 

Another possible channel to probe the pre-overlap st age 
of reionization i s studying secondary CMB anisotropics ( Hu 

199< 



supernova-d riven winds ( Madau, Ferrara, fc Rees 2001 
et al. 2001 ) or ionizations by ene rgetic X-rays ( pli 2001 



Venkatesan, Giroux, & Shull 2001). This kind of "smolder- 
ing" reionization (which results in a partial reduction in the 
neutral hydrogen fraction) will be virtually impossible to de- 
tect with radio observations. Such a signature also cannot 
be detected by measuring just the spectrum of t he flu ctua- 
tions, because, as was shown in Gnedin & Jaffe (2001), the 
power spectrum of secondary anisotropies is only weakly de- 
pende nt on the redshift of reionization. However, Gnedin & 
Jaffe (2001) noticed that the fluctuations themselves were 
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highly non-gaussian. In this paper we show how tests of the 
non-gaussianity of the secondary CMB anisotropies can de- 
tect the signature of an epoch of "smoldering reionization" . 
The physical reason behind such a possibility is simple: 
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the earlier the epoch of "smoldering reionization" begins, 
the more nonlinear the objects that are responsible for the 
produc tion of ionizing photons should be, and, therefore, the 



are known for each region in the set. Minkowski function 
als of the excursion set often called the global Mi nkowski 



functionals are particularly popular in cosmology [ 5ott et 
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highly nonlinear at early times, when the linear fluctuations 



al. 2001 



are smaller. 

A big advantage of a morphological approach is that it 
is virtually independent of assumptions about the underly- 
ing cosmological model. For example, the redshifted 21cm 
emission and absorption will depend on the rate at which 
ionization front expand into the low density IGM, gobbling 
up the neutral gas on the way. And because secondary CMB 
anisotr opics are dominated by nonlinear structures on small 
scales (Gnedin & Jaffe 2001), they are insensitive to the 
specific details of how ionization fronts expand. 



2 METHOD 
2.1 Simulations 

In this paper we use the synthetic map s of t he secondary 
CMB anisotropies from Gnedin & Jaffe (2001). These maps 
were obtained from a numerical si mulat ion of cosmological 
reionization, published in Gnedin ( |200C| ). 

The simulation included all the physical ingredients re- 
quired to study the process of reionization, including the 
3D radiative transfer. It assumed a currently fashionable 
CDM+A cosmological model with f2 m ,o = 0.3, Q.a,o = 0.7, 
h = 0.7, Qb,o = 0.04, n = 1, as = 0.91 in a comoving box 
with the size of 4/i _1 Mpc. However, for the purpose of this 
paper the specific details about a cosmological model are 
unimportant. 

For the considered cosmology the angular size of a fixed 
comoving length is almost independent of redshift for z ~ 
10. For example, our 4ft -1 Mpc box subtends 2.7 arcmin at 
z — 5, 2.2 arcmin at z = 10, and 1.9 arcmin at z = 20. Due 
to computational limitations, our highest resolution images 
include-* 4 "* 2 
from I ' 



512 pixels. Thus, we are able Lo access mullipules 
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Because our goal is to measure the degree of the non- 
gaussianity of the CMB maps, we need an appropriate tool 
for such a task. One of the most powerful mathematical tools 
to study different morphologies are Minkowski functionals. 



2.2 Minkowski Functionals 

Morphology, i.e. geometry and topology, of a random field 
can be quan tified by a set of measures know as Minkowski 
functionals (Minkowski 1903). Although being known in 
mathematics since the beginning of 20th century, they were 
brought into co smology not long ago by Mecke, Buchert & 
Wagner ( 1994) despite occasional partial rediscoveries (e.g. 
Gott et al. 1990). In the simplest form Minkowski func- 
tionals of a connected region represent three numbers char- 
acterizing the shape of a region: the area A, perimeter P, 
and genus G, which can be defined as G — 1 — where 
rih is the number of holes in the region. Minkowski func- 
tionals are additive, therefore they have precise meaning for 
any set of the regions, and can be easily computed if they 
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lytically (Longuet-Higgins 1957) 
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where R — \[2aja\ is the scale of the field; a and ai are 
the rms of the field and its first derivatives (in statistically 
isotropic fields both derivatives du/dx and du/dy have equal 
rms). It is assumed that < u >= 0. Here A(u) is the fraction 
of the total area in the regions of the excursion set, P(u) is 
the total length of the contour or total perimeter per char- 
acteristic area R 2 , and G(u) is the (number of regions) - 
(number of holes) per R 2 . 

The popularity of the global Minkowski functionals is 
probably related to the existence of the analytic formu- 
lae (eq. jjj) and computational easiness. Computing global 
Minkowski functionals is usually done simply by identify- 
ing the pixels satisfying t he threshold condition and apply- 
ing the Crofton formulae (Crofton 1968). Unfortunately, this 



simple technique is quite crude resulting in nonconvergence 
of the contour l ength. For large maps a stati stical correction 
can be applied ( Winitzki fe Kosowsky 1997 ), however its ef- 
fects in finite size maps is not quantified. The least effect 
would be a significant growth of noise. We use a different 
technique which identifies every region in the excursion set 
and computes its area, perimeter and genus from the con- 
tour points obtained by the interpolation b etween pixels. 
One can find the details of the technique in ( [Shandarin et 



al. 2002). In addition to significant increase of the accu- 



racy the method allows to use additional information about 
individual regions of the excursion set. In this study along 
with global Minkowski functionals we also use the number of 
regions (or - using the jargon of the cluster analysis - num- 
ber of clusters N c i) and the Minkowski functionals of the 
largest area region (A P ,P P and G p ) labeled as percolating 
cluster since it plays a crucial role in identifying the perco- 
lation transition. The Minkowski functionals of the perco- 
lating cluster of gaussian fields are not not known in ana- 
l ytic fo rm but can be easily computed (see e.g. Shandarin 
( 2002)). We parameterize the measurements by the total 



area of the map. The adva ntages of this parameterization 
are described in Shandarin ( | 2002| ). 

Since Minkowski functionals depend on the power spec- 
trum (R and a are obviously determined by the power spec- 
trum) of the statistical field, we need to separate this de- 
pendence from the dependence on the non-gaussian statis- 
tics. For this purpose for each simulated field of CMB 
anisotropies we produce 100 gaussian realizations with ex- 
actly the same power spectrum, and analyze these gaussian 
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Figure 1. Five Minkowski functionals and the number of re- 
gions in the excursion set as a function of fractional area for the 
simulated data (black solid lines) and 100 gaussian realizations 
(light grey dotted lines) for a 256 2 image. Medium grey long and 
short dashed lines show the median and 95% probability interval 
for the set of gaussian curves. 



realization in exactly the same manner as the simulated field. 
This also reduces the discreteness and boundary effects. We 
can then compare Minkowski functionals for the simulated 
image and for the gaussian realizations in order to measure 
the non-gaussian signal. 

An example of such a comparison for all 6 functionals 
considered here is given in Figure [[]. The behavior of all 
functionals is typical for all grids. 

Plotting the deviation from the median of 100 gaussian 
realizations better illustrates the type of non-gaussianity in 
the AT/T fields (Figure ^). Generally there are fewer re- 
gions N c i and the global perimeter P g is shorter than in the 
corresponding gaussian fields. The global g cnus Gg at very 
small A g (i.e. for rare peaks) is close to N c i, since rare peaks 
have no holes. However, at the other extreme it is signifi- 
cantly affected by the presence of holes and is different from 
the gaussian genus. The area of the percolating region does 
not differ much from the gaussian fields but the perimeter 
of the percolating region generally somewhat shorter. The 
genus of the percolating region G v shows more significant de- 
viation from the gaussian field indicating that the number 
of holes in the largest region is smaller than in the gaus- 
sian case. Comparing it with the corresponding part of the 
global genus (at A g > 0.5) one can see the influence of small 
regions: they are present in G g but omitted in G p . 

One can notice that the number of clusters Nci, 
the isocontour perimeter P g and the genus G p show the 
strongest deviations from the gaussian distribution. The 
global perimeter P g is measured in the grid units and then 
normalized to one site. 

In the rest of this paper we will focus on G p and P g 
functionals for brevity. 



Figure 2. The deviations of the simulated data (black solid lines) 
from the median of 100 gaussian realizations (light grey dotted 
lines) for a 512 2 image. Medium grey short dashed lines show the 
95% probability interval for the set of gaussian curves. 



2.3 Measuring Non-gaussianity 

The task of quantifying the degree of deviation of the simu- 
lation from the gaussian distribution is a formidable one and 
is generally an unsolved mathematical problem. Specifically, 
we would like to know how much gaussian or non-gaussian 
are the black curves in Fig. Because we cannot solve this 
problem in a general case, we have designed the following 
simple approach that allows us to measure the degree of 
non-gaussianity in some cases and identify cases where this 
cannot be done. 

We start by assuming that we can treat each curve y(x) 
in the interval < x < 1 as a collection of random variables 
yj — y(xj) where Xj,j = l...n is a subdivision of the interval 
[0, 1] (it does not have to be uniform). We then treat yj from 
gaussian realizations as gaussian random variables (this is 
the main assumption) , so that the joined probability for all 
yj for a given curve y{x) is 



i=0 j=0 

where yj = (yj) is the mean for yj and 



- Vi) 



{(.Vi - Vi)(Vi -Vi)) 



is the correlation function. Both y~j and dj can be computed 
directly from 100 gaussian realizations. 

Then for a given collection of yj from the simulated map 
we can compute the "effective" chi-square 
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The problem with this approach is that often dj is ill- 
defined if the curves are smooth and neighbouring points 
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y-j and yj+i are highly correlated. In order to avoid small 
eigenvalues in matrix Cy we use the Singular Value Decom- 
position to expand the correlation matrix as 



fc=0 



UikWkVjk, 



where Wk are eigenvalues of Cy and matrices Un, and Vjk are 
orthonormal. We can then rewrite the effective chi-square 
(0) as a sum over inverse eigenvalues, 



2 

Xoff 



n n n 

EEE 



kUjk(yi - yi)(yj -yj)- (3) 

i=0 j=0 k=0 

Let us now consider not the full sum but a partial sum, 
where we assume that Wk are sorted in the descending order 
(they are all positive since the correlation matrix is always 
positive definite), 



n n l 

Xcs(0 = ^2 EE —VikUjkiyi - m){y. 
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i=0 j=0 k=0 



The expression q(l) = Xcff(')/' as a function of / reaches 
a maximum at some value of Z max . If we interpret q(l) as 
chi-square per degree of freedom I, then Z max gives us the 
number of degrees of freedom present in the simulated curve 
yj, and we choose to truncate the sum (^) at this value 
of I - thus, discarding small u>k which makes the inverse 
correlation matrix ill-defined. 

Because g(Z max ) is not real chi-square, we cannot use 
standard tables to compute the confidence level on our sim- 
ulated curve. Instead, we compute g(Z ma x) for 100 gaussian 
curves yj and compare the cumulative distribution P(> q) 
for 100 gaussian realization with the value q for the sim- 
ulation. If q falls within the range of values that gaussian 
realizations span, we can directly estimate the probability 
that the simulation is non-gaussian by counting how many 
gaussian realizations have g(Z ma x) greater then q. However, 
because we only have 100 gaussian realization, we cannot 
measure probabilities better than 1% in this way. So, if q 
does not he within the range spanned by gaussian realiza- 
tions, we assume that g(i ma x) for gaussian realizations in- 
deed obey the chi-square distribution and compute the con- 
fidence level on q from that assumption. We also compare 
the distribution of g(7 max ) with the true chi-square distribu- 
tion using the Kolmogorov-Smirnov test, and when this test 
fails our assumption becomes invalid. In that case we can 
only say that the simulation is non-gaussian at better than 
99% confidence level. 
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Figure 3. The deviation of the genus of the percolating region 
from the median value of 100 gaussian realizations as a function 
of the total area of the excursion set for four resolutions of the 
simulated image (64 2 - 2 arcsec pixel, 128 2 - 1 arcsec pixel, 256 
- 1/2 arcsec pixel, and 512 2 - 1/4 arcsec pixel). Line markings are 
as in Fig. [l| 




Figure 4. The probability that the simulation image is a gaus- 
sian realization as a function of resolution (i.e. the image size). 
Higher resolution images are more non-gaussian. 



3 RESULTS 

Figure ^ shows the effect of angular resolution of the mea- 
surement of the Minkowski functionals. For brevity, we only 
show the genus of the percolating region as a function of 
total area of the excursion set for four values of angular 
resolution (which also implies different pixel sizes because 
the image size is fixed to 2.2 arcmin). In order to provide a 
better illustration we plot the deviation of the genus from 
the median value obtained from 100 gaussian realizations. 
As one can see, higher resolution gives progressively more 
non-gaussian signal. 



This is illustrated in Figure ^j, which shows the proba- 
bility that the signal is gaussian as a function of the image 
size. The number of clusters statistic N c i and global perime- 
ter P g show similar trends. 

This trend is easy to understand: as the resolution in- 
creases, we probe progressively smaller spatial scales that 
are more non-linear - and, thus, more non-gaussian - than 
larger scales. 

A similar trend is observed when we change the redshift 
of reionization ( we foll ow the specific procedure described in 
Gnedin & Jaffe ( [2001 )). The redshift of reionization is de- 
fined as the moment at which the rate of change of the mean 
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Figure 5. (a) deviation from the median value of the genus of 
the percolating region and (b) deviation from the median value of 
the perimeter of the excursion set as a function of area for 512 2 
resolution and four different epoch of reionization: z = 5, 7, 12, 
and 19. Line markings are as in Fig. |lj 



free path to ionizing radiation peaks. It roughly corresponds 
to the moment when the volume weighted mean fraction of 
neutral hydrogen is about 10%. 

Figure |H| shows the genus and perimeter, but other 
Minkowski functionals follow a similar trend. In this figure 
we again show the deviations from the median. 

For this result to hold, however, it does not matter 
whether the universe is reionized completely or the neutral 
fraction remains of the order of 10% or so. After all, the 
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Figure 6. The deviations of the simulated data (black solid 
lines) from the median of 100 gaussian realizations (light grey 
dotted lines) for a 512 2 image with noise. S/N =2 Medium grey 
short dashed lines show the 95% probability interval for the set 
of gaussian curves. 



CMB anisotropies are sensitive to the total optical depth. So, 
if the mean neutral fraction drops to 10% at, say, z ~ 20, and 
the complete reionization (during which the neutral frac- 
tion drops to 10~ 4 — 10~ 5 ) takes place at z ~ 6, our results 
for z — 20 would apply, whereas neither radio observations 
(because the neutral fraction is only 10%) nor optical/IR 
observations (because the neutral fraction is still 10%) will 
be able to obtain any information about the "smoldering 
reionization" period from z ~ 20 to z ~ 6. 

Real observations are almost always suffer from noise. 
In order to test the robustness of the morphological statistic 
used here, we added gaussian noise to the maps of the model 
with the reionization epoch at z = 12. Besides the model 
without noise we analyzed the models with the signal to 
noise ratios S/N = 1/5,1/4,1/3,1/2,1,2,3.3,5,10,20,50. 
We found that increasing noise from zero to about S/N = 2 
the deviations of the Minkowski functionals from the median 
of the gaussian ensemble steadily change from that shown 
in fig.^ to fig.^. The further decrease of the signal to noise 
ratio between S/N — 1/2 and S/N =1/3 eventually makes 
the non-gaussian signal undetectable as illustrated by fig.^. 
Quantitatively the probabilities of the number of regions 
N c i shown in fig.Q to be gaussian are much less than 1% 
for S/N = 2, 1 and 1/2 and is about 5% for S/N = 1/3. 
We remind that making comparison with only one hundred 
gaussian fields one cannot reliably estimate this probability 
much better then 1%. 

This conclusion is also rather insensitive to the precise 
value of the neutral fraction during the "smoldering" phase. 
Really, for our values for the cosmological parameters, the 
Thompson optical depth accumulated from the redshift in- 
terval from Z\ to Z2 is 
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Figure 7. The deviations from the median value of the number 
of regions for a 512 2 image with different ratios of the signal to 
noise. Notations are as in fig ^ 



t(zi, Zi) = cctt I n e dt 



= 2xl0- 3 [(l + ^) 3/2 -(l + ^i) 3/2 ] 



Thus, if we assume that the universe is fully reionized at 
z — 6, the interval from z = 6 to z = produces r(0, 6) = 
0.035. A period of "smoldering" reionization from z — 20 
to z = 6 with the average free electron fraction s e would 
give another contribution to the Thompson optical depth of 
t(6, 20) = 0.16a; e . For example a model in which the universe 
reionizes fully at z — 12 is very similar to the model which 
reionizes fully at z = 6 but has a period of "smoldering 
reionization" from z — 20 to z — 6 with x e = 0.35. 

Thus, until x e drops below about 0.035/0.16 = 0.2, it 
will be possible to detect the period of "smoldering" reion- 
ization. In that case, however, radio observations will likely 
be capable of tracing the neutral hydrogen abundance. 

Of course, if the period of "smoldering" reionization is 
shorter, the sensitivity to the average free electron fraction 
will be respectively weaker. How ever, for the scenario pro- 
posed in Venkatesan et al. (2001), the effect is indeed mea- 
surable. 



4 CONCLUSIONS 

We showed how the morphological analysis of the maps of 
the secondary CMB anisotropics on sub-arcminute angular 
scales can detect an extended period of "smoldering" reion- 
ization, during which the universe remains partially ionized. 
If the neutral hydrogen fraction during such a period is be- 
low about 50% but still well above 10 -5 , such a period will 
be detectable neither by radio observations of the redshifted 
21cm line nor by IR observations of the Lyman-alpha forest. 



In that case morphology of the CMB anisotropics offers the 
best chance to probe the IGM at early times. 

We computed each of six parameters as a function of 
the fractional area of the excursion set, A g : 1) the number 
of regions in the excursion set, N c i, 2) total perimeter in the 
excursion set, P g , 3) genus of the excursion set defined as the 
number of regions minus the number of holes, G g , 4) area of 
the largest (i.e. percolating) region, A p , 5) perimeter of the 
largest region, P p , and 6) the genus of the largest region, G p . 
Three parameters ( A g , P g and G g ) are known as the global 
Minkowski functionals of the excursion set, and A p , P p and 
G p are the Minkowski functionals of largest (by area) region. 
We found N c i, P g , and G p are particularly sensitive to the 
non-gaussianity in AT /T maps due to secondary reioniza- 
tion, they are also the most robust to effects of noise. A p 
and P p are not sensitive to this type of non-gaussianity, G g 
is significantly less sensitive then N c i, P g , and G p . Using N c i, 
P g , and G p one can detect the non-gaussianity in the CMB 
maps with S/N = 1/2 at the significance level of better than 
99%. 

Morphological analysis of the shapes of individual re- 
gions in the excursion set provides considerably more infor- 
mation about non-gaussianity of the maps and potentially 
may improve the current result both in terms of the resolu- 
tion of the maps and signal to noise ratio. We reserve this 
study for the future work. 
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